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Abstract 

We explore short-distance static correlation functions in the infinite XXZ chain 
using previously derived formulae which represent the correlation functions in factorized 
form. We compute two-point functions ranging over 2, 3 and 4 lattice sites as functions 
of the temperature and the magnetic field in the massive regime A > 1, extending our 
previous results to the full parameter plane of the antiferromagnetic chain (A > — 1 and 
arbitrary field h) . The factorized formulae are numerically efficient and allow for taking 
the isotropic limit (A = 1) and the Ising limit (A = oo). At the critical field separating 
the fully polarized phase from the Neel phase, the Ising chain possesses exponentially 
many ground states. The residual entropy is lifted by quantum fluctuations for large 
but finite A inducing unexpected crossover phenomena in the correlations. 

1 Introduction 

All observable information about a many-body quantum system is encoded in its correlation 
functions. For interacting systems in the thermodynamic limit this information is, in 
general, only accessible through various approximations Until rather recently very few 
examples were known, where correlation functions could be calculated exactly. Among 
the known examples there was none with short-range interactions of finite and tunable 
strength. In particular, for the large and prototypical class of integrable models that are 
solvable by Bethe ansatz, the solution was limited to the spectral properties, but except 
for what can be concluded from the finite size corrections to the spectrum rather little was 
known about their correlation functions. 
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This situation is about to change. Extensive studies by several groups, mainly focusing 
on the XXZ quantum spin chain, have eventually uncovered a structure which we believe 
to be typical for the Bethe ansatz solvable models: the one-point and the neighbour cor- 
relation functions determine everything. We call this phenomenon factorization. Based on 
the factorization we can calculate for the first time short-range correlation functions with 
arbitrary accuracy, both, for finite temperatures in the thermodynamic limit [3], and for 
arbitrary finite length in the ground state [11]. This has applications e.g. in the determina- 
tion of ESR-line shifts [28] and is useful for estimating the errors associated with standard 
numerical techniques like the quantum Monte Carlo algorithm [27]. 

Factorization appeared first [9] as the explicit factorization of certain multiple integrals 
[19,20,24] describing the density matrix elements of the isotropic Heisenberg chain in the 
ground state with no magnetic field applied. A similar factorization for the ground state 
correlation functions of the XXZ chain in the massive regime, i.e. for the model considered 
in this article, was carried out in [22,30]. After a multiple integral representation for 
finite temperature correlation functions had become available [14,15], the factorization of 
the integrals for short distances was extended to this case [5]. Following the success with 
the direct factorization of the integrals a deep exploration into the algebraic structure of 
the static correlation functions of the XXZ chain was conducted by Boos, Jimbo, Miwa, 
Smirnov and Takeyama, culminating in the works [4,7,8,21], where the factorization of all 
static correlation functions of the XXZ chain was proved in a very general setting, including 
the finite temperature and finite magnetic field cases. 

In all the recent works cited above it was crucial not to deal directly with the spin 
chain Hamiltonian, but with the associated six- vertex model [2]. The six vertex model 
naturally carries inhomogeneity parameters in horizontal and vertical directions and can 
be distorted by a disorder field without loosing its distinctive feature, the integrability. In 
this setting the density matrix of a finite segment of the XXZ chain is naturally generalized 
to include the inhomogeneity parameters and the strength a of the disorder field. The 
vertical inhomogeneity parameters Vj and the disorder field regularize the expression for the 
density matrix. The horizontal spectral parameters allow one to introduce the temperature 
into the model [14,25,26] and to adjust the Hamiltonian [31]. The density matrix of the 
inhomogeneous six-vertex model with disorder field depends polynomially on only two 
complex functions t^{y\Q) and Lo{y\^V'^a) which are efficiently described in terms of the 
solutions of integral equations [4] . The calculation of the coefficients of the polynomials is a 
purely algebraic problem. It is related to the construction of a special fermionic basis [7, 8] 
on the space of local operators acting on the space of states of the spin chain. We call this 
the algebraic part of the problem and the calculation of the functions (/9 and lo the physical 
part of the problem, since all dependence on the physical parameters, like temperature 
magnetic field or boundary conditions, is in these two functions. 

For the time being no efficient algorithm for the algebraic part of the problem is known. 
This limits the range of the correlation functions we are actually able to calculate to a few 
lattice sites. In [3] we calculated the coefficients for the two-point functions for up to 
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four lattice site by brute force computer algebra. Then we solved the physical part of the 
problem in the massless regime. At a late stage of the calculation the homogeneous limit 
and the limit of vanishing disorder field had to be taken in order to obtain the correlation 
functions of the spin chain. Due to the singular nature of this limit the first derivative of 
uj with respect to a and various derivatives with respect to the inhomogeneity parameters 
appear. The final expressions are polynomials in three functions oo, uj' (basically the a 
derivative of uj) and ip and their derivatives. The coefficients in these polynomials are 
rational functions in the deformation parameter which determines the anisotropy of the 
XXZ chain, and are algebraic and universal. Therefore the expressions derived in [3] for 
the two-point functions in the massless regime also apply in the massive regime, if one 
replaces the physical part, uj' and ip appropriately. 

This is the main subject of this work. We reformulate the expressions for the functions 
UJ, uj' and ip originally obtained in [6] in a way that is convenient for numerical calculations 
in the massive regime A > 1 of the XXZ chain. Combining this with the general formulae 
and the numerical results from [3] we obtain accurate data for the two-point functions in 
the full parameter plane of the infinite antiferromagnetic chain (anisotropy A > — 1 and 
magnetic field h arbitrary). They exhibit a surprisingly rich non-monotonic behaviour. 

We should comment on the level of mathematical rigour of our results. In [3, 6] we 
conjectured the physical part of the problem as well as the contribution to the algebraic part 
connected with the one-point functions in the limit, when a — > 0. Meanwhile we know [4,21] 
the exact physical part even for finite a and we could show [4] that it reproduces our 
conjecture for a ^ 0. The algebraic part for distances up to four lattice sites was checked 
in [3] by several independent means (comparison with the high temperature expansion of 
the multiple integrals [15], direct numerical computation, consideration of various limits). 
Now, as we may infer from [4,8,21], it is clear that the exponential formula of [3] is valid 
for arbitrary distances and zero magnetic field and at least for the distance of two lattice 
sites at finite magnetic field. On the other hand, the work [21] offers a future way for the 
rigorous calculation of the algebraic part which needs no exponential form of the density 
matrix and will hopefully lead to the calculation of static correlation functions for points 
at larger distances. 

2 Hamiltonian and density matrix 

We consider the Hamiltonian 

N 

j=-N+l 

of the XXZ chain in the massive antiferromagnetic regime (J > and A > 1). The 
cr", J = —N + 1, . . . ,N, act locally as Pauli matrices, A = ch(?/) is the anisotropy para- 
meter and J the exchange coupling. 
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The XXZ Hamiltonian preserves the z-component of the total spin 

TV 



Thus, the magnetization in z-direction is a thermodynamic quantity, and the thermal 
equilibrium of the finite system is characterized by the statistical operator 

g-(3<jv-^Sjv)/T 

PN(T,h) = TTTT ,o -./rp (3) 

depending on the temperature T and the external magnetic field h. 

In the thermodynamic limit L = 2N oo the system severely simplifies, since, from 
the six-vertex model point of view, a single eigenstate of the so-called quantum transfer 
matrix determines the state of thermodynamic equilibrium and hence all static correlation 
functions [14] . Clearly the naive limit N ^ oo makes no sense in ([3]) . To perform this limit 
in a sensible way we fix integers m, n with m < n and introduce the density matrix 

D\m,n\{T,h) = lim tr_7v+l,...,m-l,n+l,n+2,...,Ar /5Ar(r, /l) (4) 

N—i-oo 

of the segment [m, n] of the infinite chain which is well defined. It has the reduction 
properties 

t^m D[jn,n] (T, h) = D[m+l,n] ^) i -^[m,n] {T, h) = -D[m,n-1] {T, h) . (5) 

We consider the vector space W of operators on the infinite chain which act non-trivially 
only on a finite number of lattice sites (c|<7| for arbitrary fixed j and k is an example of 

such an operator). On this space we define a map trfr,;„„i 

trf [m,n] = ■ • • ^ tr.m-2 \ trm_i \ tr„+i \ tr.„+2 . . . . (6) 

It restricts the action of to the interval [m,n]. 
Because of ^ the definition Z : W — > C, 

Z(0) = lim trm,...,n (7", ^)trfr„ „i (7) 

<r<~i — ^ — rv^ L ' J L ' J 



-OO 

n^oo 



makes sense and determines the thermal expectation value of 0. The functional Z has the 
natural interpretation of the statistical operator on the space W and may be thought of as 
the thermodynamic limit oi pN, equation 

The functional Z allows for a generalization within the framework of the six-vertex 
model with disorder field. In the seminal paper [21] this generalization was denoted Z'^ 
(with K, referring to the magnetic field in certain units). It is this functional which depends 
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on only two functions ip{i'\a) and U2\oi) and which clearly separates the physical and 
the algebraic part of the problem. For details we refer the reader to [4,21]. 

In this work we are dealing with applications. We shall only need an understanding of 
the general idea of factorization and the concrete formula obtained in [3]. For this reason 
we refrain in the following from any sophisticated mathematical notation and simply write 
(0)/i,T = Z{0) for thermal expectation values on the infinite chain. 

3 The physical part of the construction 

As explained in [3] the limit a ^ leads to three functions uj{fii,fi2), ^' ifJ-i, and ip{fJ,) 
that determine the correlation functions of the XXZ chain. A relatively simple description 
of these functions in terms of an auxiliary function a and a generalized magnetization 
density G was given in [6]. We call it the a- formulation. The a- formulation is useful for 
deriving multiple integral formulae [14] and for studying the high-temperature expansion of 
the free energy and the correlation functions [32] . It allows for a rather uniform description 
of the massless and the massive XXZ chain, the only difference between the two cases being 
a different choice for the so-called canonical integration contour in the complex plane [14]. 

For numerical calculations of thermodynamic properties [25,26] or short-distance corre- 
lators from the multiple integral [3,10] one has to change to a different formulation we refer 
to as the bb-formulation. It needs pairs of functions, but the defining integral equations 
involve only straight contours. It is here where one has to distinguish between the massless 
and the massive case, when it comes to numerical calculations. Two separate computer 
programs are needed. In [3] we studied the massless case. Here we proceed to the massive 
case. We start from the o-formulation of the physical part proposed in [6] and switch to 
the bb-formulation. This is a standard procedure which (in the massive case) basically re- 
quires the application of Fourier series and the convolution theorem. As far as the integral 
equations are concerned the reader may e.g. refer to [10]. Then it is rather obvious how 
to proceed in a similar way with the expressions for the functions tp, uj and iv' (equations 
(15), (17) and (20) for a = in [6]). For this reason we just state the results. 

Let us define a pair of auxiliary functions b , b as the solution of the non-linear integral 
equations (NLIE) 

In b(x) = - ^ - '^^^-P^d{x) + K * In (1 + b) (x) - * In (l + b) (x) , (8a) 

In b(x) =^ - ^^^-p^d{x) + K * In (1 + b) (x) - K+ * In (1 + b) (x) (8b) 

with / * g{x) = ^ '^yfi-^ ~ y)9{y) denoting a convolution. 

The integral equations are specified by the integration kernels k, and the driving 
terms, which contain a single transcendental function d. Note that the physical parameters 
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temperature T, coupling J and magnetic field h enter into the calculation of the correlation 
functions only through the driving terms of the NLIE ([SD . 

The functions d and well as some other functions occurring in other integral 

equations below, have simple Fourier series representations which arise naturally in the 
derivation of the non-linear equations and which are also useful for solving them numeri- 
cally, 



d{x)= 



°° ^i2kx 



ch{r]k) 



°° ^-ri\k\+2ikx 



2 ch.(r]k) 



K^{x) = k{x lb irj^) . 



(9a) 
(9b) 



Alternatively the kernel and the function d can be realized in terms of special functions. 
For the function d we find 

2K , /2Kx .7]^ 



d{x) 



dn 



\ vr vr 



(10) 



where dn(x,r) is one of the Jacobi elliptic functions [1], and K is the complete elliptic 
integral of the first kind. 

The integration kernel k can be expressed in terms of the normalized trigonometric 
gamma function T introduced in [29] . In order to define it we first of all recall the definition 
of the g-gamma function [13], 



T,{z) = (1 - qf-^ n r 



1-g" 

_ „z+n+l 



n=l 



Then, for q = exp(— 4//), 



T{r-z) 



r,(i/2) 

Tq{iz + 1/2) 



exp 



1 , rz^ 

- — mvr H izm 

2 2 



2r 



and 



nix) = — dx In 



r(2^;|^)r(2r?;-^-|) 



2r, 



r(2,?;-|-)r(2r,;|--§) 



(11) 



(12) 



(13) 



In addition to the auxiliary functions b and b we need two more pairs of functions 
and g' ^ to express cu, u>' and 93 by means of auxiliary functions. Both are solutions of linear 
integral equations involving b and b. 



K * 



9,1 
1 + b- 



1 + b- 



j{x). 



g^ (x) = -d{x - n) + K* Y^^ix) - K+ * Y^ppi (a;) 



(14a) 

(14b) 
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and 



a'ti^) = - VC+{x - fi)+r]l* Y^^^ix) -r]l * ^ (x) 



9't . . - 



+ * Y^(^) - * iri^^^) • (15b) 

The new integration kernel / and the new functions occurring in (jlSh are conveniently 
described as Fourier series, 

fc=— oo ^ ' ' 

where we used the standard convention sign(O) = 0, and 

°° ±r]k+2ikx 
fe=— oo ^ ' ' 

The functions a;(//i, /Lt2), a;'(/ii,/i2) and which determine the inhomogeneous cor- 
relation functions can be expressed as integrals over the above functions. The function 

dx ( 9~^{x) glix) \ 

— 7r/2 

witli jl = —ifi is related to the magnetization 

m(/i,r) = i(a|>^^^ = -1(^(0) (19) 
which is the only independent one-point function of the XXZ chain. The function 

u;{^l,^i2) = -M^^2-ill) + K^{i:i2-ili)-d* + (;i2) (20) 

with 

K (x) = (21) 

2 sin(x -|- irj) sm{x — irj) 



^In the following this convention will also be used for fii,fL2. 
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also determines the internal energy [3]. The last function which is necessary to determine 
the correlation functions is 

-^^-*TT|^(/^2)-r?c+*-i^(/i2), (22) 



where 



W = , v- (23) 

2 sm(x + zr/j sin(x — i?7) 

The physical meaning of is less intuitive. Its derivative appears in the two-point functions 
and is therefore related to certain neighbour correlators. 

As shown in [3,6] it is necessary to perform the homogeneous limit for the calculation of 
correlation functions. This limit seems rather singular, because the coefficients coming from 
the algebraic part in general have poles, when two inhomogeneities coincide. Still, these 
poles are canceled by zeros coming from certain symmetric combinations of the functions 
uj{fii, 112), i^'(/^i,/^2) and fin) in the numerators. In order to perform the limit one has to 
apply I'Hospital's rule. This finally leads to polynomials in the functions a;(0, 0), uj'{0,0), 
<p(0) and in derivatives of these functions with respect to the inhomogeneity parameters 
evaluated at zero. We denote derivatives with respect to the first argument by subscripts 
X and derivatives with respect to the second argument by subscripts y and leave out zero 
arguments for simplicity. Then e.g. coxyy = dxdyLo{x,y)\x^y=o etc. 

For the examples in the next section the non-linear integral equations for b and b as 
well as their linear counterparts for gj^^ and g'^^^^ were solved iteratively in Fourier space, 

using the fast Fourier transformation algorithm. The derivatives of g^^ and g'^^^^ with 
respect to fi, needed in the computation of the respective derivatives of to and to' , satisfy 
linear integral equations as well, which were obtained as derivatives of the equations for 
gj^^ and g'^^^ ■ The modifications due to the derivatives are particularly simple in Fourier 
space. 



4 Examples of short distance correlators 

In the sequel we shall restrict ourselves to the longitudinal and transversal two-point func- 
tions {o'fa^)j, ^, c^) J- /j for n = 2, 3, 4H The algebraic part for these correlation functions 

■^Our method is based on the density matrix and aUows us to obtain aU its matrix elements. In particular, 
we may calculate all other independent two-point functions like e.g. (ffcrJOj,^ for n — 2,3,4. Still, since 
this will give more bulky expressions of the type ([26}, p7p . we refrained from the temptation of presenting 
them all. 
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was already calculated in our previous papers [6] and [3]. Here we merely cite those results. 
For the nearest neighbour correlation functions the following equations were obtained 



Kf^l)r,h=cth(r/)u; + ^ 



T,h 



2sh(77) 

The expansions for next-nearest neighbours read 



2ri 



{^"i^Dtm =2cth(2r/)w + + 



uj'x , t'h{r]){uJxx - 2uJxy) sh'^i'nWxxy 



I „X X 



1 



sh(2??) 



ch(2??) , ch(2?/)th(r/)(u;^^ - 2a;a;j^) ^\?{ri)J^. 



2i] 



-UJ, 



+ 



■'xxy 



8r] 



(24a) 
(24b) 

(25a) 
(25b) 



The length of the formulae grows rapidly with the number of lattice sites. For n = 4 it 
reads 



1 



|384g4 (1 + (5 _ 4^2 ^ 5^4^ ^2^ 

- 8 (1 + (52 + 64g2 _ 234^'' + 64q^ + 52q^ + q^'^)) if 



UJ. 



xy 



+ 192g^ 



-I + q^f {I + Aq^ + q^) ifujyy 



+ (-1 + ^2)^(1 + /) (1 + 4^2^/) ,^2 



IQSq^ {-I - q^ + q^ + q^) i]J y 



+ 16 (-1 + q^f (1 + 6g2 + Ug^ + n^e ^ g^8 ^ ^10^ 



xyy 



2{-l + q^f{l + 2q^ + 2q' + q') r^co',. 



2 ' 2q^ 

,4 



xxyyy 



LOy — UJUJ. 



+ 8 (-1 + q^f (1 + g2) (1 + + Uq^ + Gg^ + g^) rf 
+ (-1 - 4g2 - 22g^ - 12g'5 + 12g^° + 22^^^ ^ 4^14 ^ ^16^ ^2 

-j- \2UJyy(jJxy ~t~ AujyLOyyy 1 2u^y AoJOJxyyy ~t~ GUJU^xXyy 



xy 
2 



-6w, 



10 



UJyyUJ y LOyUJ yy "t" UJUJ 



xyy 



+ (-1 + q^f (1 + + ^ ^ ^ 



^^xyyy^ y ^^xxyy^ y '^^yyy^ yy 



-\- QUJrj^yyUJ yy "t" lUJ yyUJ yyy r^yUJ yyy QUJyyUJ ^yy "t" AUJyUJ Xyyy '2UJUJ XXyyy 



+ 3(-l+(?^)'(l+g2 + g4) 

and in the transversal case, 

1 



+ ijj' ^Uj' 



xxyyy 



T,h 



3072^5 (-1 + g6)r^2 

- 768g^ (1 + 10g2 + 7?^^ 



+ 16g2 (-1 + g^)^ (31 + 56^^ - 30g^ + 56g^ + 31.?^) rj^- 



96g2 (-1 + q^f (3 + hq^ - Aq^ + hq^ + 'iq^) r] 



8\ J2 



yy 



^^xyyy ^'^^xxyy 



+ q2 (_l + g2)4^^^4^2^^4) ^2 

+ 192g2 (-3 -q^-q^ + q^ + glO + Sr?^^) rjLo'y 
+ 8 (-1 + g^)^ (1 - 12^2 _ 25/ - 25^*5 - 12q^ + 7]uj[ 
+ 2{-l + q^)'{l + 2q^ + 2q^ + q^) tju;' 



xyy 



xxyyy 



+ 16g2 (-1 + g2)3 ^ ^ ^ ^7^6^ 



.6\ „2 



^xy — ^y 



+ q^ (-5 - V - ISg"^ + 13g^ + 4g^° + 5g^2) r?^ 



12a;yj^ - 2AL0yyL0xy 



8uJyU!yyy + 24^^ (jJajj^y + SoJUJxyyy — i2uiUJxxyy 



-l + g2)4(l_9/ 



+ (-1 + g^)^ (1 + + ^ 5^6 ^ ^8^ ^ 



iOyyUJ y — UJyLO yy + WW ^J/J/ 

-AlOxyyy^ y ~1~ ^^xxvv^ ' 



^xxyy^ y 



(26) 



-|- ^.LOyyyLJ yy 6CJ^y|yCJ |y|; 'ZtOyyiO yyy AuJ xijU) yyy "h 6U^|y|;U^ x^/ly AiOyLO XVVV ~^ '2^UJUJ XX' 



^xyy^ yy 



+ 3(-l + /)'(l + <?2 + /) 



^yyU/ yyy 



xy^ yyy 



^yy^ xyy 



oj yyy^ xyy ~\~ ^ yy^ ^yyy 



^y^ xyyy 



y^ xxyyy 



xxyyy 

, (27) 



with q = . 

Using the representations via Hnear and non-Unear integral equations for the functions 
Lo and Lo' of the previous section we can determine high-precision numerical values for the 
various two-point correlators. Figures [T]S] show selected examples of connected two-point 
functions. The longitudinal correlation functions (o"f(T^)y^ — {(^i) t h i'^n) t h (^ ~ 2,3,4) 
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are shown in the left panels, while the right panels show the transversal correlation functions 
(crfcj^)^^ — {^n)Th- Note that the contributions from the one-point functions are 

somewhat trivial. Due to the translational invariance of the Hamiltonian the longitudinal 
one-point functions do not depend on the site index, {(tDj.^ = (cr^)^^;^, and are given by 
(jlSp for all n. The conservation of the z-component of the total spin, on the other hand, 
implies that {crn)Th ~ ^- "^^^ above notation for the connected two-point functions is 
merely used for systematic reasons. 

In figure [T] we show the dependence of the correlation functions on the magnetic field 
for A = 2 and for several values of the temperature. At low temperatures one can clearly 
see the two critical fields of the ground state phase diagram, figure [5l The saturation field 
is at /i = 12 and the critical field at which the excitation gap opens can be calculated 
[12, 33] to be located at /i ~ 1.55921. This can be seen even better in figure [H where 
the {crfal)rp^ — (crf)^,^ WDrh correlation is shown for different low temperatures and for 
magnetic fields in the vicinity of the two critical fields. Here the curves for T/J = 0.01 
and T/J = 0.001 are clearly distinguishable, which is not the case on the larger scale of 
figure [TJ For fields above the saturation field h = 12 the correlators have a fixed value as 
saturation sets in. For h below the lower critical field we observe a constant behaviour as 
well, which in this case is due to the excitation gap. In between, however, for values of the 
magnetic field for which the system is critical, an interesting non-monotonic behaviour can 
be seen. It is most pronounced for the correlators (o"fo"|)^^ — {(tDj, ^ {^Dt h which two 
local extrema exist if the temperature is sufficiently small. If the temperature is too high, 
thermal fiuctuations dominate and all correlations die out. In figure [1] this is illustrated 
with the curve for T/J = W in the n = 2 case. For n = 3, 4 the effect is similar. 

In figure [2] we show the correlation functions for several values of the magnetic field 
as functions of the temperature. Here we observe a non-monotonic behaviour as well 
for intermediate magnetic fields. Regarding the low temperature behaviour we notice 
that the connected correlation functions depend in a non-trivial way on the distance and 
on the magnetic field. There exists a field for which the connected correlation functions 
have maximal modulus. For nearest neighbours the figures for the transversal correlation 
functions on the right side, show the largest values for h/J = 6, but for next nearest and 
next-to-next nearest neighbours the correlations for h/J = 10 are more pronounced at low 
temperature. For the field strength h/J = 16, which is above the saturation field, the 
connected correlation functions deviate the lesser from zero at intermediate temperatures 
the larger the distance is. 

The excitation gap can also clearly be seen in the left panel in figure [3l where we show 
i^i^Dxh ~ ^'^Dt h i'^Dxh temperatures and several magnetic fields in the vicinity 

of the lower critical field. With increasing magnetic field the correlations start deviating 
from the h = line at lower and lower temperatures, until finally the zero temperature 
limit is different from the zero field case. Close to the critical field the system is particularly 
sensitive to small changes in the field strength, which can be well observed by comparing 
e.g. the correlation functions for h = 1.56, slightly above the critical field, and h = 1.55, 
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slightly below the critical field. 

The right panel again shows (cfcrDj,^ — (o'f')j,^ i'^Dxh^ ^^^^ time as a function of 
A for fixed small temperature T/J = 0.01 and several magnetic fields. For large values 
of A the connected correlation function is independent of the magnetic field due to the 
excitation gap. The effect of the phase transition from the massive fully polarized state 
to the critical antiferromagnet is also visible in this figure. For h/J = 8, for instance, we 
observe an abrupt change at A = 1, which, at this field strength, is the critical anisotropy 
at which saturation sets in. Note that for producing the data for < A < 1 in the figure we 
used the formulation and computer implementation of our previous work [3]. The curves 
smoothly match at A = 1. 
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Figure 1: (Color online) (crftr^) — {o'i){a^) and {crfa^) — (o"f)(o"J^) for different values of 
temperature T/J and fixed anisotropy A = 2. The rows are for n = 2, 3, 4. 
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Figure 2: (Color online) (crftr^) — {o'i){a^) and {crfa^) — (o"f)((T^) for different values of 
the magnetic field h/J and fixed anisotropy A = 2. The rows are for n = 2,3, 4. 
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Figure 3: (Color online) The left panel shows (crfal) — {crl){crl) for different magnetic fields 
h/J and fixed anisotropy A = 2 at low temperatures. The right panel shows {afal) — 
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Figure 4: (Color online) (crfcr|) — (cif) ((t|) for different low temperatures T/J and fixed 
anisotropy A = 2. The values of the magnetic field are in the vicinity of the critical fields. 
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Figure 5: (Color online) Ground state phase diagram of the XXZ chain. 
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5 Ising limit 

The subject of this work is the study of correlation functions of the XXZ chain for A > 1. 
This includes the two extreme cases A = 1 and A = oo. The physically interesting 
isotropic point A = 1 is difficult to access analytically. In fact, the regularization of the 
density matrix by a disorder field mentioned in the introduction fails exactly at this point, 
and it seems that the picture of factorization described above needs a slight modification in 
the presence of a magnetic field, when additional independent functions, called 'moments' 
in [5], appear. On the other hand, the numerics discussed in the previous section remains 
remarkably stable, if one approaches the isotropic point from above, and smoothly matches 
the numerics for approaching it from below [3]. 

For A = oo or, more precisely, in the limit A ^ oo for finite J/ = J A the XXZ 
Hamiltonian ([1]) turns into the Hamiltonian of the Ising chain 



N 



"Ki = Jj ^ 



h 



N 



(28) 



-N+1 



-N+1 



The Ising chain may be viewed as a classical one-dimensional model of statistical mechanics 
with 2 X 2-transfer matrix (see ()39p below). For this reason its correlation functions can 
be calculated explicitly [2], 



sh (2^] 



sh2(A) +e4>/^/^ 



sh^ (A) 



sh" 



'A 

'.2T 



+ 



ch (22,) 



sh2(A) +e4J,/T 



) + ^'''^^ ^^^2 (A) + e^Ji/T ^ /sh2(^)+eWT 



(29a) 
(29b) 



An interesting and non-trivial test of the formulae of the previous section is to reproduce 
these results analytically and numerically in the limit A — > 00. A similar exercise starting 
from the multiple integral formulae of [14] was carried out in [17]. In our case at hand we 
start with the auxiliary functions in the bb-formulation which turn into 



+ Wsh" 



h 

2r 



+ e 



4J//T 



b{x) 



^h/2T 



-sh I 



2T) 



+ Jsh^ (A) + e^^^/T 



(30a) 
(30b) 
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in the Ising limit. Note that, in this special case, they do not depend on the spectral 
parameter. This is also true for the functions 



-1 ± 



sh(A) 



sh^ (A^ 



which moreover do not depend on fi anymore. The same is therefore true for 

sh(A) 



(31) 



and 



'^(Aii,At2) 



cth 



sh2(A) +e4j,/r 



2Ji 



+ 



ch(A) 



sh 



sh2 (A) + e^Ji/T 



(32) 



(33) 



Similarly, for g' , for which a rescaling is necessary to obtain a sensible limit, one obtains 



1 



1 + 



ch(A) 



V 



/sh2(A) +e4J//r^ 

Note that this rescaling is unproblematic for 77 — > oo due to (fTSj) . Thus, we obtain 

w'(//i,//2) 



(34) 



0. 



V 



(35) 



The above leads to the correct results for (cjf) and (crfcrl)- For the next-to-nearest 
neighbour and all higher correlation functions, however, some of the coefficients in the 
factorized form of the correlation functions may diverge in the Ising limit, e.g. the term 
sla^ {rj)/r] on the right hand side of (I25ap . All the terms where such type of divergence 
occurs have to be handled separately. 

Sticking with the example of equation ()25aP we define 



/ (x) = lim sh{rj)d g- (x) 

?7— >oo ^ 

h^{x) = lim sh(r/)9^^ ^ 



?7— >oo 



11=0 



(36a) 
(36b) 
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implying 

f^i^)=-Y^e''''-^e-'''' (37a) 

/-(a;)=-e*2.__^g-i2x (37b) 

h+i^)=-~e^''' (37c) 

h-ix)=j^e-^'^ (37d) 

with b, b according to (j30p . Using the latter expressions we conclude that 



^-i;-^, ^ ^^y (i + a)(l + b) 



hm^ — = 4 - , , (38) 



which indeed reproduces (|29bp when inserted into the equation (j25ap for (afo"!). 

After having shown that the Ising limit is included in our representation of the correla- 
tion functions by means of solutions of NLIE, it appears natural to compare the analytic 
limit with numerical results for large A. In figure [6] we show the temperature dependence of 
the connected two-point function (afa^) — {af) ((t|) for the Heisenberg chain near the Ising 
limit, for different values of the anisotropy parameter and the magnetic field. In general, 
these curves do not depend on A for large anisotropics and match the curves for the Ising 
limit given by (j29p . E.g., the curves in the left panel for A = 1000 and h/Jj = 3.9 or 
h/Jj = 4.1, respectively, match the corresponding Ising curves to the precision of the line 
width in the figure. In the vicinity of /i = 47/ , however, we observe large deviations from 
the low-temperature Ising curves. 

This behaviour is induced by a critical point at 1/A = 0, h = 4 J/ in the zero temper- 
ature phase diagram of the XXZ chain (see figure [7|) which separates regions of different 
asymptotics for the two-point correlation functions of the Ising chain. The latter fact is 
most easily understood by inspection of the 2x2 transfer matrix of the Ising chain [2] , 

g/i/2T g2J^/T 



* - g2J//T g-V2T J • (39) 

Its eigenvalues X± (A+ > A_) and eigenvectors v± determine the partition function of the 
Ising chain as well as its two-point functions. In fact, when we wrote (f29bl) . we used the 
formula [2] 

i^K+i) = (v+,ctV+)2 + (v+,a^v„)2 (40) 

and the explicit expression for the eigenvalues and eigenvectors of t, which the reader may 
easily calculate from (|39p . In order to understand the different zero temperature behaviour 
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of the Ising correlation functions it suffices to consider the low temperature asymptotics 
of the transfer matrix. Using ([5U|) , ([Klj) one easily distinguishes three different asymptotic 
regimes. 

. h>AJr. o)> (af<+i>~l, (41a) 

. h<AJi: i ~ e^^^/^ ( ? ) ' ^ ' ^^^^^ 

Thus, in the Ising limit the absolute value of the ground state correlation functions 
is generally independent of the spatial separation, except at the critical field he = 4:Jj, 
where we see exponential decay due to the residual entropy resulting from an exponential 
degeneracy of the ground state. For finite A this degeneracy is lifted by residual quantum 
mechanical interactions causing algebraically decaying correlations. The critical point of 
the Ising chain corresponding to a first order phase transition appears as a critical point 
in the ground state phase diagram of the XXZ chain if we draw it in the h-l/A plane for 
fixed J[ , see figure [71 

Contrary to the Ising case, there are always two critical fields if A is finite. They 
correspond to two second order phase transition lines at the boundary of the critical phase 
in the h-l/A diagram. This means that, away from the Ising limit, even the ground state 
correlation functions depend continuously on the magnetic field. Hence, for sufficiently 
low temperatures, there must be deviations from the Ising curves for values of h and A 
belonging to the critical phase. In general the zero temperature limit should depend on 
the anisotropy and the magnetic field. This can be seen in the left panel in figure [6] for 
h/Jj = 3.9965, A = 1000 and A = 2000. However, as the width of the critical phase gets 
smaller with increasing A, one can find for each magnetic field, except for the critical field, 
a sufficiently large A such that for larger anisotropics the correlation functions behave as 
in the Ising model. 

For h/Jj = 4, there is always a deviation in the zero temperature behaviour between the 
Ising chain and the XXZ chain for finite A. See figure El where the connected correlation 
functions are shown for A = 1000, 2000, 5000 and in the Ising limit. Interestingly in 
this case the zero temperature limit does not depend on the (finite) value of A. This is a 
general behaviour as we observe that the correlation functions are asymptotically constant 
on straight lines ending in the Ising critical point. We can see this exemplarily in the right 
panel of figure [6l where the connected correlation functions show the same zero temperature 
asymptotics for three different pairs of values of the anisotropy and the magnetic field. 

Comparing these curves with those for the Ising model with the same magnetic fields, 
we identify four different temperature regimes. For very high temperatures the curves are 
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independent of A and h. Then, for intermediate temperatures, they are independent of 
the (large) anisotropy and only depend on the magnetic field. This is where the curves 
for finite A and the curves of the Ising model match. Next comes an regime where the 
correlation functions depend on the anisotropy and the magnetic field, and finally, for very 
low temperatures, only the product {h — hc)A determines the correlation functions. In 
fact, this is a rare example, where the full low temperature asymptotics of thermodynamic 
quantities and short-range correlation functions in the vicinity of a critical point can be 
worked out analytically [16]. 
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Figure 6: (Color online) (cjf (t|) — (af ) (cr|) for large values of the anisotropy, fixed Jj and 
different values of h/Jj. The labels in both panels are the tupels A,h/Jj where 'Ising' 
denotes the analytic Ising curves. 
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Figure 7: (Color online) Ground state phase diagram of the XXZ chain in the /i-l/A plane 
for fixed J/ = J A and large A. 
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6 Paramagnet 

Another limit, which is interesting from a technical point of view, is the paramagnet J = 0. 
At a first glance it might look rather trivial. Yet, as the coefficients in section 2] depend on 
7], the functions to, to' and ip depend on r] as well, but the final result must not. In this sense 
the paramagnetic limit is even more intricate than the Ising limit, as, in contrast to the 
latter, all the summands in the algebraic part contribute. This way we have an additional 
test for the coefficients in (IM1) - (|27|) . We note that the physical part takes the form 



in the paramagnetic limit. 

7 Conclusion 

We have studied the short-range correlation functions of the XXZ chain in the massive 
phase by means of the equations ()24p - (j27p representing them in factorized form. For this 
purpose we derived in section [3] a representation of the physical part of the correlation 
functions which is well suited for the implementation on a computer. We obtained high- 
accuracy data for the correlation functions of the spin chain in the thermodynamic limit. 
Together with our previous results [3] we can now access the full parameter plane of the 
infinite antiferromagnetic chain (anisotropy A > — 1 and magnetic field h arbitrary) at 
arbitrary temperatures. The short-range correlation functions show a surprisingly rich non- 
monotonous behaviour at intermediate magnetic fields and temperatures. At low enough 
temperatures the critical lines of the ground state phase diagram can be read off from our 
data. The numerics is stable in the isotropic limit A 1+ and in the Ising limit A — > co. 

So far our method is still limited in the accessible range of the correlation functions. In 
order to extend this range, an efficient algorithm for the calculation of the algebraic part 
of the correlation functions is needed. Since this is a well defined, purely algebraic problem 
though, we have little doubt that it will be solved in the future. For the exact calculation 
of the large-distance asymptotics at finite temperatures, on the other hand, we believe that 
new insight will be needed (for recent progress on the corresponding ground state problem 




(42b) 



(42a) 



(42c) 



see [23]). 
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